/*
 heat-spectral lib. solves time-dependent PDEs using SEM-NI method.
 Copyright (C) 2009, Altieri Michele and Stefanoni Andrea
 Copyright (C) 2008, Proverbio Alessandro and Radice David

 This program is free software: you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
 the Free Software Foundation, either version 3 of the License, or
 (at your option) any later version.
 
 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 GNU General Public License for more details.
 
 You should have received a copy of the GNU General Public License
 along with this program. If not, see <http://www.gnu.org/licenses/>.
*/

/**
   \file integrator.cpp
   \brief Gaussian integration class

    This class implements the Gaussian integration on the GLL nodes on a given
    rectangular dominium. An affine transformation is used to transoform the
    domain in the square \f$ [-1,1] \times [-1,1] \f$.
*/

# include <cmath>
# include <cassert>

# include "integrator.hpp"

namespace spectral {

Integrator::Integrator(int nnodes,
                       double delta,
                       double eps,
                       int nmax,
                       int order)   : _weightsX(NULL),
                                      _weightsY(NULL),
                                      _nnodes(nnodes),
                                      _nodes(NULL),
                                      _weights(NULL),
                                      _delta(delta),
                                      _eps(eps),
                                      _nmax(nmax),
                                      _order(order) {
    // Using less then 4 nodes would require another constructor
    // or a slower one. Neither the first nor the latter is desired.
    assert(_nnodes>3);

    _weightsX=new double[_nnodes];
    _weightsY=new double[_nnodes];

    _nodes=new double[_nnodes];
    _weights=new double[_nnodes];
    Poly L(_nnodes-1);
    legendre(L);
    Poly dL(L.differentiate());
    double * lnodes=horner(dL,_delta,_eps,_nmax);

    // Sorting the nodes
    std::qsort(lnodes,_nnodes-2,sizeof(double),&(_compare));

    _nodes[0]=-1;
    _nodes[_nnodes-1]=1;
    for(int i=1;i<_nnodes-1;++i) {
        _nodes[i]=lnodes[i-1];
    }
    for(int i=0;i<_nnodes;++i) {
        _weights[i]=2.0/double(_nnodes*(_nnodes-1));
        assert(std::abs(L(_nodes[i]))>1e-16);
        _weights[i]=_weights[i]/std::pow(L(_nodes[i]),2);
    }
    delete[] lnodes;

    // ==============================================
    // Lagrange Polynomials Generation
    // ==============================================

    Poly tmp;
    int n=_nnodes-1;

    // Initializate Polynomials
    _lagrange=new Poly * [_order+1];

    // ----------------------------------------------
    // First Order Polynomials
    // ----------------------------------------------
    _lagrange[0]=new Poly[n+1];
    for(int i=0;i<n+1;++i){
        _lagrange[0][i].setdegree(n);
    }

    // First Polynomial
    double alpha=(n*(n+1))*L(-1.0);
    _lagrange[0][0][0]=-dL[0]/alpha;
    _lagrange[0][0][n]=dL[n-1]/alpha;
    for(int i=1;i<n;++i) {
        _lagrange[0][0][i]=(-dL[i]+dL[i-1])/alpha;
    }
    // Middle Polynomials
    for(int i=1;i<n;++i) {
        alpha=(n*(n+1))*L(_nodes[i]);
        tmp=deflate(dL,_nodes[i]);
        _lagrange[0][i][n]=tmp[n-2]/alpha;
        _lagrange[0][i][n-1]=tmp[n-3]/alpha;
        _lagrange[0][i][0]=-tmp[0]/alpha;
        _lagrange[0][i][1]=-tmp[1]/alpha;
        for(int j=2;j<n-1;++j) {
            _lagrange[0][i][j]=(-tmp[j]+tmp[j-2])/alpha;
        }
    }
    // Final Polynomial
    alpha=(n*(n+1))*L(1.0);
    _lagrange[0][n][0]=dL[0]/alpha;
    _lagrange[0][n][n]=dL[n-1]/alpha;
    for(int i=1;i<n;++i) {
        _lagrange[0][n][i]=(dL[i]+dL[i-1])/alpha;
    }
    // ----------------------------------------------

    // ----------------------------------------------
    // Derivatives polynomials
    // ----------------------------------------------
    for(int k=1;k<=_order;++k) {
        _lagrange[k]=new Poly[_nnodes];
        for(int i=0;i<_nnodes;++i) {
            _lagrange[k][i]=_lagrange[k-1][i].differentiate();
        }
    }
    // ----------------------------------------------

    // ==============================================

}

Integrator::Integrator(Integrator const & integr) : _nnodes(integr._nnodes),
                                                    _delta(integr._delta),
                                                    _eps(integr._eps),
                                                    _nmax(integr._nmax),
                                                    _order(integr._order)
{
    _weightsX=new double[_nnodes];
    _weightsY=new double[_nnodes];

    _nodes=new double[_nnodes];
    _weights=new double[_nnodes];

    _lagrange=new Poly * [_order+1];
    for(int i=0;i<=_order;++i) {
        _lagrange[i]=new Poly[_nnodes];
    }

    for(int i=0;i<_nnodes;++i) {
        _nodes[i]=integr._nodes[i];
        _weights[i]=integr._weights[i];
    }

    for(int i=0;i<=_order;++i) {
        for(int j=0;j<_nnodes;++j) {
            _lagrange[i][j]=integr._lagrange[i][j];
        }
    }
}

Integrator::~Integrator() {
    delete[] _weightsX;
    delete[] _weightsY;
    delete[] _nodes;
    delete[] _weights;
    for(int i=0;i<=_order;++i) {
        delete[] _lagrange[i];
    }
    delete[] _lagrange;
}

double Integrator::integrate(double const * v,
                             Point<double> const & A,
                             Point<double> const & B) const {
    double lengthX=B.get(0)-A.get(0);
    assert(lengthX>0);
    double lengthY=B.get(1)-A.get(1);
    assert(lengthY>0);

    // Precalculating weights and nodes is more efficient
    // it's O(n) instead of O(n^2):
    for(int i=0;i<_nnodes;++i) {
       _weightsX[i]=0.5*lengthX*_weights[i];
       _weightsY[i]=0.5*lengthY*_weights[i];
    }

    double sum(0.0);

    for(int i=0;i<_nnodes;++i) {
        for(int j=0;j<_nnodes;++j) {
            sum+=_weightsX[i]*_weightsY[j]*v[i+j*_nnodes];
        }
    }

    return sum;
}

TensorProd Integrator::lagrange(int i) const {
	assert( i>=0 && i<_nnodes );
	Poly ones(1.);

	return TensorProd(_lagrange[0][i],
			  ones,
			  i,
			  0,
			  0,
			  0,
			  1.0,
			  1.0);
}

TensorProd Integrator::lagrange(int i, int j) const {
    assert( i>=0 && i<_nnodes );
    assert( j>=0 && j<_nnodes );

    return TensorProd(_lagrange[0][i],
                      _lagrange[0][j],
                      i,
                      j,
                      0,
                      0,
                      1.0,
                      1.0);
}

TensorProd Integrator::Dx(TensorProd const & tp, double lengthX) const {
    double weightX=tp._weightX;
    double weightY=tp._weightY;
    int i=tp._idx;
    int j=tp._idy;
    int odx=tp._odx+1;
    int ody=tp._ody;

    assert( odx <= _order );
    assert( ody <= _order );

    return TensorProd(_lagrange[odx][i],
                      _lagrange[ody][j],
                      i,
                      j,
                      odx,
                      ody,
                      2.0/lengthX*weightX,
                      weightY);
}

TensorProd Integrator::Dy(TensorProd const & tp, double lengthY) const {
    double weightX=tp._weightX;
    double weightY=tp._weightY;
    int i=tp._idx;
    int j=tp._idy;
    int odx=tp._odx;
    int ody=tp._ody+1;

    assert( odx <= _order );
    assert( ody <= _order );

    return TensorProd(_lagrange[odx][i],
                      _lagrange[ody][j],
                      i,
                      j,
                      odx,
                      ody,
                      weightX,
                      2.0/lengthY*weightY);
}

int Integrator::getnumnodes() const {
    return _nnodes;
}

double Integrator::getnode(int n) const{
    assert( n>=0 && n<_nnodes );
    return _nodes[n];
}

double * Integrator::getnodes() const {
    double * nodes;
    nodes=new double[_nnodes];
    for(int i=0;i<_nnodes;++i) {
        nodes[i]=_nodes[i];
    }
    return nodes;
}

int Integrator::_compare(const void * a, const void * b) {
    double x( *(double*)a - *(double*)b );
    return (x<0)?(-1):(1);
}

} // namespace
